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Abstract 

Hybrid  renewable  energy  systems  are  expected  to  become  competitive  to  conventional  power  generation  systems  in  the  near  future  and,  thus, 
optimization  of  their  operation  is  of  particular  interest.  In  this  work,  a  hybrid  power  generation  system  is  studied  consisting  of  the  following  main 
components:  photovoltaic  array  (PV),  electrolyser,  metal  hydride  tanks,  and  proton  exchange  membrane  fuel  cells  (PEMFC).  The  key  advantage  of 
the  hybrid  system  compared  to  stand-alone  photovoltaic  systems  is  that  it  can  store  efficiently  solar  energy  by  transforming  it  to  hydrogen,  which 
is  the  fuel  supplied  to  the  fuel  cell.  However,  decision  making  regarding  the  operation  of  this  system  is  a  rather  complicated  task.  A  complete 
framework  is  proposed  for  managing  such  systems  that  is  based  on  a  rolling  time  horizon  philosophy. 

©  2007  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Over  the  last  few  years,  the  limited  reserves  in  conventional 
fuels  like  oil  and  the  increasing  interest  of  the  public  in  the 
protection  of  the  environment,  have  led  the  research  commu¬ 
nity  to  seek  for  alternative  Renewable  Energy  Sources  (RES). 
Special  emphasis  has  been  given  on  the  development  and  imple¬ 
mentation  of  fuel-cell  systems,  for  both  academic  purposes  and 
industrial  applications.  Fuel  cells  may  be  considered  as  continu¬ 
ous  chemical  reactors  which  convert  fuel  and  oxidant  chemical 
potential  into  electrical  energy.  The  key  advantages  of  fuel  cells 
compared  to  the  conventional  electrical  power  generation  tech¬ 
nologies  are:  higher  efficiency,  especially  when  the  waste  heat 
is  used  for  co-generation,  quiet  operation  suitable  for  residen¬ 
tial  applications,  and  almost  zero  levels  of  produced  pollutant 
gases.  These  advantages  are  due  to  the  fact  that  power  genera¬ 
tion  in  fuel  cell  systems  is  not  based  on  combustion  techniques 
and  temperature  gradients.  Thus,  the  limitations  imposed  by  the 
second  law  of  thermodynamics  on  the  operation  of  fuel  cells 
are  much  less  severe  than  the  limitations  imposed  on  conven¬ 
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tional  energy  conversion  systems,  which  are  traditionally  used 
for  power  generation.  The  most  important  drawback  concerning 
fuel  cell  technology  is  that  the  basic  fuel  they  use  (hydrogen) 
does  not  exist  free  in  nature  [1]. 

A  mature  technology  of  electrical  power  generation  from 
RES  is  the  photovoltaic  (PV)  technology,  which  is  character¬ 
ized  by  zero  emissions  of  pollutant  products  but,  also,  by  low 
efficiency  levels.  However,  low  efficiency  should  not  be  consid¬ 
ered  as  a  major  disadvantage,  due  to  the  zero  cost  of  solar  energy. 
The  main  disadvantage  of  the  PV  technology  is  that  availability 
of  solar  energy  is  limited  to  daylight  periods  and  depends  on  the 
geographical  location  and  the  weather  conditions.  It  is  obvious 
that  solar  energy  cannot  be  stored,  unless  it  is  transformed  to 
another  type  of  energy.  Furthermore,  the  accurate  prediction  of 
solar  energy  distribution  is  not  feasible  over  a  long  future  time 
horizon  [2]. 

In  order  to  overcome  the  aforementioned  limitations,  the  con¬ 
cept  of  combining  RES  with  hydrogen  technologies  emerged; 
and,  systematic  studies  of  this  innovative  type  of  combined  sys¬ 
tems  appeared  in  the  early  90  [3].  The  term  “hybrid  power 
generation  system”  refers  to  all  systems  that  combine  differ¬ 
ent  energy  technologies  (RES,  Hydrogen,  Biomass  etc.)  in 
order  to  meet  the  required  electrical  and  thermal  loads  of  the 
consumer. 
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Nomenclature 

Ai  anisotropy  index 

Cost  cost  for  purchasing  electricity  from  the  grid 

(€  kWh-1) 

CostCool  cost  per  unit  volume  of  H2  for  removing  heat 
from  the  hydride  tanks  (€  Nm-3) 

CostHeat  cost  per  unit  volume  of  H2  for  adding  heat  to 


the  hydride  tanks  (€  Nm“3) 

cs 

energy  demand  (kW) 

eff 

inverter  efficiency  (%) 

EFG 

electricity  from  the  grid  (kW) 

EFIn 

electrolyzer  input  power  (kW) 

EFOut 

hydrogen  production  rate  from  the  electrolyzer 
(Nm3  h-1) 

ETG 

electricity  to  the  grid  (kW) 

/ 

factor  used  to  account  for  horizon  brightening 

FCIn 

fuel  cell  hydrogen  consumption  (Nm3  h-1) 

FCOut 

fuel  cell  power  generation  (kW) 

GSIT 

global  solar  irradiance  on  tilted  surface 

4 

beam  irradiance  (kW  m-2) 

4 

diffuse  irradiance  (kW  m-2) 

Inv 

hydrogen  stored  in  hydride  tanks  (Nm3) 

INVIn 

inverter  input  power  (kW) 

INVOut 

inverter  output  power  (kW) 

GSI 

global  solar  irradiance  (kW  m“2) 

GSIT 

global  solar  irradiance  on  tilted  surface  (kW  m-2) 

Price 

profit  from  selling  electricity  to  the  grid 
(€  kWh-1) 

PDir 

direct  transferred  power  (kW) 

PVOut 

photovoltaic  output  power  (kW) 

P  V surface  net  surface  area  of  the  photovoltaic  array  (m  2 ) 

Rb 

ratio  of  the  beam  irradiance  on  the  tilted  surface 
of  the  photovoltaic  array  to  the  beam  irradiance 
on  the  horizontal  surface 

4mb 

ambient  temperature  (°C) 

Greek  symbols 

p 

slope  of  the  photovoltaic  array  (°) 

Pg 

ground  reflectance 

Abbreviations 

AGGSI 

adjusted  Gaussian  global  solar  irradiance 

FC 

fuel  cells 

GSI 

global  solar  irradiance 

MINFP 

mixed  integer  non  linear  programming 

MPC 

model  predictive  control 

NNM 

neural  network  model 

PEM 

proton  exchange  membrane 

PV 

photovoltaic 

RBF 

radial  basis  function 

RES 

renewable  energy  sources 

RESHS 

renewable  energy  system  with  hydrogen  storage 

Many  researchers  have  focused  their  studies  on  such  systems 
by  following  three  main  directions:  (a)  experimental  studies 
on  established  systems;  (b)  simulation  and  optimization  of 
their  operation  and  (c)  real  time  control.  A  number  of  stud¬ 
ies  [4-7]  present  the  experimental  experience  which  was  gained 
during  the  operation  of  hybrid  systems  (solar-hydrogen,  wind- 
hydrogen)  in  pilot  scale.  Representative  references  along  the 
second  direction,  on  the  other  hand,  are  also  available  [3,8-17]. 
Santarelli  and  Macagno  [8],  analyzed  a  solar-hydrogen  system 
for  a  residential  building  in  the  Alps.  The  purpose  of  that  system 
was  to  satisfy  the  complete  demand  in  electrical  energy  and  part 
of  the  consumer  heat  demand.  Kelouwani  et  al.  [9]  developed  a 
dynamic  model  for  a  stand-alone  wind/solar-hydrogen  system 
whereas  Korpaas  et  al.  [10]  proposed  a  method  for  scheduling 
and  operation  of  wind  power  plants  in  electricity  markets.  Indica¬ 
tive  studies  concerning  the  simulation  and  the  optimization  (in 
terms  of  choosing  the  optimal  main  components  for  establishing 
the  whole  system)  of  a  renewable  energy  system  with  hydro¬ 
gen  storage  (RESHS)  are  given  in  references  [11-17].  Finally, 
certain  references  [18-20]  present  different  methodologies  of 
controlling  hybrid  renewable  energy  power  generation  systems 
for  given  load  and  meteorological  conditions. 

In  this  study,  a  novel  framework  is  proposed  for  the  real-time 
operation  optimization  of  a  solar-hydrogen  power  generation 
system.  The  framework  is  built  in  several  steps:  In  the  first  step, 
a  neural  network  model  (NNM)  for  predicting  global  solar  irradi- 
ance  (GSI)  distribution  on  horizontal  surfaces  is  developed.  Next 
step  is  the  estimation  of  the  electrical  energy  produced  by  the  PV 
array.  Then  a  model  is  developed  that  describes  realistically  the 
performance  and  the  limitations  of  the  hybrid  system.  Finally,  the 
rolling  horizon  fundamental  principle  of  Model  Predictive  Con¬ 
trol  (MPC)  is  used  for  developing  an  optimal  decision  strategy. 
In  particular,  an  optimization  problem  is  formulated  based  on  the 
discretization  of  a  future  prediction  horizon,  which  determines 
the  optimal  values  of  the  decision  variables  at  each  discrete  time 
instance.  The  optimization  problem  is  reformulated  and  solved 
at  each  particular  time  instance,  thus  allowing  the  model  to  be 
updated  with  new  information  as  the  process  progresses  in  time. 

The  rest  of  the  paper  is  structured  as  follows:  The  next  section 
describes  the  hybrid  energy  system  that  is  studied  in  this  work. 
Section  3  presents  step-by-step  the  development  of  the  decision 
making  framework,  from  the  generation  of  the  GSI  prediction 
model  to  the  formulation  of  the  complete  optimization  prob¬ 
lem  which  provides  the  optimal  values  of  the  decision  variables. 
Section  4  presents  simulated  results  of  the  proposed  system  for 
a  typical  household  and  illustrates  how  the  cost  for  purchasing 
electrical  energy  can  be  reduced.  Finally,  the  paper  ends  with 
the  concluding  remarks  of  Section  5. 

2.  System  description 

A  detailed  schematic  of  the  particular  hybrid  system  is  shown 
in  Fig.  1.  The  notation  used  in  this  diagram  and  throughout  the 
rest  of  the  paper  is  explained  in  the  Nomenclature. 

The  main  components  of  the  RESHS  are  the  Photovoltaic 
Array  (PV),  the  Electrolyzer,  the  Metal  Hydride  Tanks,  the  Pro¬ 
ton  Exchange  Membrane  Fuel  Cells  (PEMFC)  and  the  inverter. 


PL.  Zervas  et  al.  /  Journal  of  Power  Sources  181  (2008)  327-338 


329 


Photovoltaic 

Array 


from  the  grid 


Fig.  1.  Configuration  of  the  solar-hydrogen  system. 


The  system  operates  as  follows:  the  electrical  energy  generated 
from  the  PV  may  either  supply  the  electrolyzer  -  in  order  to  pro¬ 
duce  hydrogen — or  cover  directly  the  consumer  requirements. 
The  hydrogen  produced  is  stored  in  metal  hydride  tanks  con¬ 
nected  to  a  fuel  cell  stack.  Operation  of  the  fuel  cells  produces 
electrical  energy,  which  is  used  to  satisfy  part  or  all  of  the  con¬ 
sumer  needs.  Additionally,  the  system  gives  the  opportunity  to 
purchase  electricity  from  the  grid  or  sell  electricity  to  the  grid. 

3.  Development  of  the  optimization  procedure 

In  the  first  step,  a  NNM  for  predicting  the  GSI  distribution 
on  horizontal  surfaces  is  developed.  Next  step  is  the  estimation 
of  the  electrical  energy  produced  by  the  PV  array.  Then  a  model 
is  developed  which  realistically  describes  the  performance  and 
the  limitations  of  the  hybrid  system.  Finally,  the  rolling  horizon 
fundamental  principle  of  MPC  is  used  for  designing  an  optimal 
decision  strategy. 

3.1.  Local  prediction  of  the  global  solar  irradiance 

Prediction  of  GSI  is  essential  for  the  proper  design  of  build¬ 
ings  and  solar  energy  systems.  In  particular,  as  far  as  PV  systems 
are  concerned,  estimation  of  GSI  and  of  its  components  (beam, 
diffuse,  and  albedo  irradiance)  is  necessary  for  sizing  purposes 
[2] .  In  the  optimization  framework  described  in  this  paper,  pre¬ 
diction  of  GSI  locally  (where  the  hybrid  system  lies)  is  crucial, 
because  it  provides  the  necessary  information  for  estimating  the 
PV  power  generation  over  a  future  horizon.  GSI  is  predicted 
using  a  novel  methodology  that  is  based  on  neural  network  tech¬ 
nology  and,  in  particular,  on  the  Radial  Basis  Function  (RBF) 
neural  network  architecture  [21]. 

The  GSI  daily  distribution  is  approximated  by  a  Gaussian- 
type  function.  The  RBF  neural  network  is  used  to  predict  the 
parameters  of  the  Gaussian-function  (amplitude  and  width), 
using  the  weather  state  and  the  daylight  duration  as  input  vari¬ 
ables.  The  weather  conditions  are  categorized  into  discrete 
states,  whereas  the  daylight  duration  is  represented  in  daylight 
tenths.  The  main  meteorological  parameter  that  is  used  to  clas¬ 
sify  the  weather  state  is  the  presence  and  the  type  of  clouds. 
In  accordance  with  many  weather-forecast  models,  six  different 
states  were  defined  for  the  classification  of  cloudiness,  namely: 


•  State  1:  Rainfall 

•  State  2:  Heavy  clouds 

•  State  3:  Cloudy 

•  State  4:  Partly  cloudy 

•  State  5:  Few  clouds 

•  State  6:  Clear 

Fig.  2  presents  characteristic  daily  solar  irradiance  distri¬ 
butions  for  all  six  weather  states.  The  NNM  is  developed  by 
applying  the  fuzzy  means  RBF  training  algorithm  [22]  on  a 
database  which  contains  local  observations  of  the  input  vari¬ 
ables  (weather  condition  and  the  daylight  duration)  as  well  as  the 
parameters  of  the  Gaussian  function  over  a  period  of  time  (one 
calendar  year  is  strongly  recommended).  A  correction  method¬ 
ology  for  the  two  tails  of  the  Gaussian  function  improves  further 
the  accuracy  of  the  produced  model.  The  prediction  model  that 
is  finally  developed  is  named  adjusted  Gaussian  Global  solar 
irradiance  model  (AGGSI).  This  mathematical  model  is  used 
to  provide  reliable  future  predictions  of  the  daily  GSI  distribu¬ 
tions  on  horizontal  surfaces,  given  only  the  weather  state  forecast 
and  the  daylight  duration.  The  procedure  for  generating  an  esti¬ 
mate  of  the  GSI  distribution  for  a  particular  day  is  presented 
graphically  in  Fig.  3  and  is  described  in  details  in  [21]. 

3.2.  Estimation  of  the  electrical  power  generation  of  the  PV 
array 

Assuming  that  the  model  described  in  the  previous  subsection 
is  available,  the  clearness  index  can  be  calculated  as  the  ratio  of 
the  GSI  to  the  extraterrestrial  irradiance.  The  latter  is  computed 
through  astronomical  and  geographical  data,  given  the  solar  con¬ 
stant,  the  day-of-year  and  the  time-of-day.  The  components  of 
GSI  on  the  horizontal  plane  (beam  irradiance,  4  and  diffuse  irra¬ 
diance,  If)  are  estimated  next,  using  a  correlation  equation  which 
expresses  the  diffuse  fraction  (ratio  of  the  diffuse  irradiance  to 
the  GSI)  as  a  function  of  the  clearness  index.  Among  several 
equations  that  are  available  in  the  literature,  the  Orgill  and  Hol¬ 
lands  [23]  correlation  was  selected  as  the  most  reliable  for  the 
region  of  Athens,  Greece.  In  order  to  estimate  Global  Solar  Irra¬ 
diance  on  Tilted  surface  (GSIT),  it  is  necessary  to  employ  three 
additional  parameters:  (a)  the  geometric  factor,  R b,  defined  as 
the  ratio  of  beam  irradiance  on  the  tilted  surface  of  the  PV  array 
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Fig.  2.  Characteristic  daily  GSI  distributions  for  all  six  weather  states. 


to  beam  irradiance  on  the  horizontal  plane;  (b)  the  anisotropy 
index,  Aj,  which  is  a  measure  of  the  atmospheric  transmittance 
of  beam  irradiance  and  is  defined  as  the  ratio  of  beam  irradiance 
to  the  extraterrestrial  irradiance,  and  (c)  the  modulating  factor, 
/,  which  is  used  to  account  for  “horizon  brightening”  and  is  cal¬ 
culated  as  the  square  root  of  the  beam  irradiance  to  GSI  ratio. 
Given  the  ground  reflectance  pg,  and  the  slope  of  the  PV  array 
ft,  GSIT  can  be  estimated  through  the  following  equation  ([2], 
§2.16): 


GSIT  =  (7b  +  /d  At)Rh  +  /d(l  -  A/) 


1  +  COS  ) 


x  (l  +  /sin3  (0)+GSIpg(l — 2—^)  (1) 


Fig.  3.  The  proposed  algorithm  for  predicting  the  daily  GSI  distribution. 


Finally,  an  empirical  model  is  developed  which  expresses 
the  photovoltaic  power  generation  (PVOut)  as  a  function  of 
GSIT  and  the  ambient  temperature  (Tamb),  using  historical  data. 
The  recorded  measurements  were  collected  during  two  years  of 
operation  of  the  Polycrystalline-Si  PV  arrays  integrated  on  the 
southern  fa$ade  and  the  roof  of  the  NTUA  Chemical  Engineer¬ 
ing  Building  [24].  The  above  data  set  was  used  to  produce  an 
empirical  model  which  correlates  the  aforementioned  quantities. 
Our  study  showed  that  the  following  linear  equation  provides  an 
excellent  correlation  for  the  array  DC  power  generation: 

PVOut  =  (0.128  GSIT  -  0.239  x  10“3  Tamb)  x  PVsurface 

(2) 

where  PVsurface  is  the  net  surface  area  (not  including  the  frames 
of  the  panels)  of  the  PV  array.  The  units  of  the  Tamb  factor  are 
in  (kW  m-2  °C_1). 

The  above  equation  was  developed  and  validated  using  a 
random  partition  of  the  available  data  into  a  training  set  and 
a  validation  set  in  a  ratio  75%:25%.  Only  the  training  data  were 
utilized  for  the  development  of  the  model  and  the  rest  of  the  data 
were  used  for  validation  purposes.  The  coefficient  of  determi¬ 
nation  R 2  (0.9909)  corresponding  to  the  validation  data  set  is 
representative  of  the  high  correlation  between  experimental  and 
predicted  values.  In  order  to  show  that  the  success  of  the  model 
was  not  the  result  of  a  chance  correlation,  100  different  models 
were  developed  based  on  100  different  random  partitions  of  the 
data  into  training  and  validation  sets.  Robustness  of  the  model 
was  clearly  illustrated,  since  the  coefficients  of  determination 
R2  were  over  0.989  for  all  100  partitions. 

In  conclusion,  PVOut  can  be  predicted  locally  over  a  future 
time  horizon  given  only  standard  weather  forecast  information 
(weather  state  and  ambient  temperature).  According  to  the  pro¬ 
duced  mathematical  model,  the  effect  of  the  ambient  temperature 
on  the  PV  performance  is  negative  (i.e  when  this  temperature 
increases,  the  efficiency  of  the  array  drops).  This  observation 
agrees  with  the  results  that  have  been  reported  in  the  litera- 
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Fig.  4.  Estimation  of  PV  power  generation. 


ture  (cf.  [2],  ch.  23).  Fig.  4  illustrates  the  procedure  used  for 
estimating  the  PV  power  generation  over  a  future  time  horizon. 

3.3.  Estimation  of  the  inverter  efficiency 

The  inverter  efficiency  is  a  function  of  the  normalized  input 
power,  i.e.  the  ratio  of  the  inverter  input  power  to  the  nominal 
input  capacity.  However,  for  normalized  input  power  greater  than 
15%,  the  efficiency  for  most  commercial  inverters  is  close  to  0.9, 
as  it  is  verified  by  the  manufacturer  specifications.  In  our  model, 
we  assume  that  the  inverter  efficiency  is  constant  over  the  entire 
range  of  normalized  input  power.  This  assumption  is  also  used  in 
the  popular  TRNSYS  simulator  [25].  In  order  to  compute  exper¬ 
imentally  an  estimate  of  the  inverter  efficiency,  we  used  again 
historical  data  that  were  collected  over  the  course  of  two  years 
operation  from  a  commercial  inverter  located  in  the  Solar  Engi¬ 
neering  Unit  of  Chemical  Engineering  School  at  NTUA  [24]. 
The  estimated  efficiency  was  found  equal  to  89.3%.  The  result 
was  validated  using  the  same  procedure  that  was  described  in 
the  previous  subsection,  i.e.  by  partitioning  the  data  into  train¬ 
ing  and  validation  sets  in  a  ratio  75%:25%.  We  generated  100 
different  random  partitions  and  in  all  of  them  the  coefficient  of 
determination  R 2  was  over  0.99. 

3.4.  Development  of  the  mass  and  electrical  energy  balance 

The  two  main  processes  of  the  hybrid  system  occur  in  the 
electrolyzer  and  the  fuel  cells,  while  the  hydride  tanks  serve 
as  storage  units.  The  inverter  efficiency  is  also  taken  into 
account  in  the  development  of  the  equations  that  are  presented 
in  this  subsection.  Using  the  notation  CS(k),  EFG(&),  ELIn(^), 
ELOut(k),  ETG(k)  FCIn(jfc),  FCOut(ifc),  INVIn(k),  INVOut(&), 
PDir(k),  PVOut(^),  to  denote  the  integrated  values  of  the  vari¬ 
ables  introduced  in  Fig.  1  over  the  length  of  a  particular  time 
period  k  and  the  notation  Inv(k)  for  the  hydrogen  inventory  level 
at  the  end  of  period  k,  we  can  write  the  relevant  equations. 

3.4.1.  Electrolyzer  and  fuel  cell  performance  equations 

Manufacturers  of  electrolyzers  provide  the  necessary  data 
to  generate  equations  that  express  the  power  consumed  by  the 
electrolyzer  as  a  function  of  the  produced  rate  of  hydrogen. 
Similarly,  for  fuel  cell  systems  equations  can  be  developed  that 


relate  the  produced  power  with  the  consumed  rate  of  hydrogen. 
These  equations  take  into  account  the  parasitic  loads,  which 
are  due  to  the  auxiliary  subunits  in  each  particular  system.  The 
equations  differ  between  manufacturers,  so  at  this  point  we  use 
generic  functions  to  denote  the  existence  of  such  mathematical 
relationships. 

ELIn(k)  =  Li(ELOut(k))  (3) 

FCOut(k)  =  L2(  FCIn(ifc))  (4) 

3.4.2.  Hydrogen  mass  balance 

The  hydrogen  inventory  level  at  the  end  of  period  k  is  equal 
to  the  inventory  level  at  the  end  of  the  previous  time  period,  plus 
the  amount  of  hydrogen  produced  by  the  electrolyzer  during 
period  k  minus  the  quantity  of  hydrogen  consumed  by  the  fuel 
cell  during  the  same  period. 

Inv(k)  =  Inv(k  -  1)  +  ELOut(k)  -  FCIn(k)  (5) 

3.4.3.  Overall  electrical  energy  balance 

The  overall  electrical  energy  balance  takes  into  account  all  the 
different  areas  in  the  system  where  energy  is  generated  or  con¬ 
sumed.  In  order  to  arrive  to  this  equation,  we  first  write  electrical 
energy  balance  equations  at  nodes  1,  2,  3,  which  are  indicated 
in  Fig.  1. 

Energy  balance  at  node  1  : 

PVOut(k)  +  EFG(k)  =  ELIn(fc)  +  PDir(k)  (6) 

Energy  balance  at  node  2  :  PDir(k)  +  FCOut(k)  =  INVIn(k) 

(7) 

Energy  balance  at  node  3  :  INVOut(k)  =  CS(k)+ETG(k)  (8) 

If  we  additionally  assume  a  constant  efficiency  inverter  eff, 
as  indicated  in  the  previous  subsection,  i.e. 

INVOut(k)  =  eff  x  INVIn(k)  (9) 

then  from  Eqs.  (6)-(9)  we  can  easily  derive  the  following  overall 
electrical  energy  balance  equation. 

eff  x  (PVOut(k)  +  EFG(k)  +  FCOut(k)  -  ELIn(k)) 

=  CS  (k)  +  ETG(k) 


(10) 


332 


PL.  Zervas  et  al.  /  Journal  of  Power  Sources  181  (2008)  327-338 


3.5.  Formulation  of  the  optimization  problem 

MPC  has  recently  advanced  to  a  popular  control  methodology 
for  industrial  and  process  applications.  Its  wide  adoption  from 
the  industry  is  largely  based  on  the  inherent  ability  of  the  method 
to  handle  efficiently  constraints  and  uncertainties  in  multivari¬ 
able  dynamical  systems.  MPC  is  based  on  the  following  simple 
idea:  at  each  discrete  time  instance  the  control  action  is  obtained 
by  solving  on-line  a  finite-horizon  open-loop  optimal  control 
problem,  using  the  current  state  of  the  system  as  the  initial  state. 
A  finite-optimal  control  sequence  is  obtained,  from  which  only 
the  first  element  is  actually  applied  to  the  system.  The  procedure 
is  repeated  after  each  state  transition. 

The  fundamental  rolling  horizon  principle  of  MPC  is  used  in 
this  work  to  design  an  optimal-decision  strategy  for  the  RESHS. 
The  model  developed  in  the  previous  subsection  is  extended  to 
account  for  the  future  time  periods,  while  the  estimated  profiles 
of  the  electrical  energy  produced  by  the  PV  array  and  the  energy 
demand  provide  the  input  information.  An  optimization  prob¬ 
lem  is  formulated  based  on  the  discretization  of  the  future  time 
horizon.  The  solution  of  the  optimization  problem  provides  the 
optimal  values  of  the  decision  variables  over  the  future  horizon. 
The  key  decision  variables  are:  the  energy  which  is  supplied  to 
the  electrolyser,  the  energy  purchased  or  sold  to  the  grid  dur¬ 
ing  each  future  time  period,  and  the  quantity  of  hydrogen  stored 
in  the  hydride  tanks  at  the  end  of  each  time  period.  The  adop¬ 
tion  of  the  rolling  horizon  principle  means  that  only  the  optimal 
values  associated  with  the  first  future  time  period  are  actually 
applied  to  the  system.  The  model  is  updated  with  new  informa¬ 
tion  as  the  process  progresses  with  time  and  the  optimization 
problem  is  reformulated  and  solved  at  each  particular  discrete 
time  instance. 

Assuming  a  number  of  P  equally  spaced  future  time  inter¬ 
vals,  the  precise  formulation  of  the  optimization  problem  is  as 
follows: 

Objective  function: 

EFG(&  +  i)  x  Cost(&  +  i)  —  ETG(&  +  i)  x  Priced  +  i) 
+Inv  U(k  +  i)  x  CostCool  +  Inv  D(k  +  i)  x  CostHeat 

(ID 

Constraints: 

If  Inv  (A:  +  i)  -  Inv(k  +  i  -  1)  >  0 

then  Inv  U(k  +  i)  =  [Inv(k  +  i)  -  Inv(k  -fi -1)1  (12) 

i  =  1, ...,  P 

If  Inv(k  +  i )  —  Inv(k  +  i  —  1)  <  0 
then  Inv  D(k  +  i)  =  [Inv(k  +  i  —  1)  —  Inv(k  +  /)],  (13) 

i  =  1, ...,  P 

ELIn(£  +  i)  =  Li(ELOut(k  +  /)),  i  =  1,  ...,  P  (14) 

FCOut(k  +  i)  =  L2( FCIn(k  +  /)),  i  =  1,  ...,  P  (15) 

Inv(k  +  i )  =  Inv(&  +  /  —  1)  +  ELOut(&  +  i)  —  FCIn(k  +  /), 
i  =  l,...,P  (16) 


eff  x  (PVOut(k  +  i)  +  EFG(k  +  i)  +  FCOut(k  +  i ) 

— ELIn(£  +  /))  =  (17) 

CS (k  +  i)  +  ETG(k  +  /),  i  =  1,  P 

If  EF G(k  +  0  >  0  then  ETG(&  +  0  =  0,  i  =  1 , ...,  P 

(18) 

If  ETG(&  +  0  >  0  then  EFG(k  +  i)  =  0, 
i  =  l,...,P  (19) 

Eq.  (11)  is  the  objective  function  to  be  minimized,  which 
actually  represents  the  cost  for  purchasing  electricity  during  the 
future  time  horizon  minus  the  corresponding  profit  of  selling 
electricity  to  the  grid.  Cost(&  +  i)  and  Price  ( k  +  i)  represent  the 
cost  per  unit  and  profit  per  unit  of  purchasing/selling  electricity 
from/to  the  grid  within  time  period  k+i.  Thus,  the  system  is 
able  to  take  into  account  situations  where  electricity  prices  vary 
from  period  to  period  (consider  for  example  the  different  day 
and  night  electricity  rates).  The  objective  function  also  takes  into 
account  the  cost  for  removing  heat  from  the  hydride  tanks  during 
addition  of  hydrogen  and  the  cost  for  adding  heat  to  the  tanks 
during  extraction.  The  cost  is  considered  proportional  to  the 
amount  of  hydrogen  ln\U(k  +  i)An\D(k  +  i)  added  to  or  extracted 
from  the  hydride  tanks.  In  particular  the  amount  of  hydrogen 
added  during  a  time  interval  is  multiplied  by  CostCool,  while  the 
amount  of  hydrogen  removed  is  multiplied  by  CostHeat,  where 
CostCool  and  CostHeat  represent  the  cost  per  unit  volume  of 
hydrogen  for  removing/adding  heat  to  the  hydride  tanks.  As  far 
as  the  constraints  are  concerned,  Eqs.  (12)  and  (13)  compute 
the  amounts  of  hydrogen  added  or  subtracted  from  the  hydride 
tanks  during  time  intervals,  which  are  used  in  the  objective  func¬ 
tion.  Eqs.  (14)— (17)  are  extensions  of  Eqs.  (3)-(5),  (10)  for  the 
entire  future  time  horizon.  Eqs.  (18)  and  (19)  prevent  the  system 
from  purchasing  and  selling  electrical  energy  to  the  grid  dur¬ 
ing  the  same  time  interval.  In  addition  to  the  above  constraints, 
the  variables  are  bounded  between  upper  and  lower  limits  to 
account  for  process  limitations.  PVOut(k  +  i),  ELin(/r  +  i),  ELOut 
(k  +  0,  FCIn(k  +  i)  and  FCOut(k  +  i)  are  bounded  from  above  by 
the  production  capacities  of  the  PV  array,  the  electrolyzer  and 
the  fuel  cell  system,  while  Inv(k  +  i)  is  bounded  from  above  by 
the  storage  capacity  of  the  hydride  tanks.  All  the  variables  are 
obviously  bounded  from  below  by  zero.  Alternatively,  we  can 
impose  greater-than-zero  lower  bounds  to  prevent  the  system 
components  from  operating  at  very  low  production  rates  and 
the  hydrogen  inventory  level  from  going  below  a  safe  threshold 
value.  Additional  constraints  may  be  included  in  the  optimiza¬ 
tion  problem,  according  to  each  particular  situation. 

Remark:  The  optimization  problem  described  by  Eqs. 

(1 1) — (19)  and  the  upper  and  lower  bounds  on  the  decision 
variables  is  a  complex  nonlinear  programming  problem.  Eqs. 

( 1 2) -(  13)  and  ( 1 8)-(  1 9)  have  been  presented  in  an  intuitive  form. 
However  insertion  of  the  if-then  statements  in  mathematical  pro¬ 
gramming  solvers  is  not  allowed  or  may  cause  computational 
flows.  The  same  equations  can  be  formulated  so  that  they  do 
not  make  use  of  if-then  statements,  by  introducing  sets  of  binary 
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Formulation  and  solution 
of  a  Mixed  Integer  Nonlinear 
Programming  Problem 


Fig.  5.  Formulation  of  the  optimization  problem  applied  to  the  RESHS. 


variables.  In  particular,  Eqs.  (12)  and  (13)  can  be  formulated  as 
follows: 

Inv(k  +  i)  —  Inv(k  +  i  —  1)  =  Inv  U(k  +  i )  —  Inv  D(k  +  /), 
i  =  1, P 

Inv  U(k  +  i)  <  Mi  x  IC(£  +  /),  i  =  1,  P 

Inv  D(k  +  i)  <  Mi  x  IH (k  +  /),  i  =  1, P 

I C(k  +  i)  +  IH (k  +  i)  =  1,  i=  1,  P 

(20) 

where  Mi  is  a  fixed  big  positive  number  (the  storage  capacity  of 
the  hydride  tanks)  and  IC (k+i),  IH (k  +  i),  1=1,. .  .,P  are  binary 
variables. 

Similarly,  Eqs.  (18)  and  (19)  can  be  written  in  the  following 
form 

EF G(k  +  i)  <  M2  x  JP(k  +  /),  i  =  1, ...,  P 

ET G(k  +  i )  <  M2  x  IS (k  +  /),  i  =  1,  ...,  P  (21) 

IP(k  +  i )  +  IS  (it  +  i)  =  1,  i  =  1, ...,  P 

where  M2  is  a  fixed  big  positive  number  (an  upper  bound  on 
the  amount  of  energy  that  may  be  purchased/sold  from/to  the 
grid  during  a  time  interval)  and  IP(&  +  i),  IS {k  +  /),/=  1,. .  .,P  are 
additional  sets  of  binary  variables. 

Using  the  formulations  described  by  Eqs.  (20)  and  (21)  in 
place  of  Eqs.  (12)— (13)  and  (18)— (19)  a  Mixed  Integer  Nonlinear 
Programming  problem  (MINLP)  is  generated.  The  problem  is 
solved  using  GAMS  [26],  which  is  a  powerful  commercial  mod¬ 
elling  language.  The  formulation  of  the  optimization  problem  is 
described  graphically  in  Fig.  5. 

4.  Results 

In  this  section,  simulation  results  are  presented,  concerning 
the  application  of  the  proposed  method  in  optimizing  the  perfor¬ 
mance  of  a  RESHS,  which  is  used  to  cover  the  energy  needs  of  a 
typical  household  located  in  Athens,  Greece.  Based  on  recorded 
data,  reference  energy  consumption  daily  profiles  were  created 


for  each  month.  Fig.  6  shows  the  reference  profile  for  October, 
which  will  be  used  for  the  simulations  that  follow  in  the  sequel. 

Future  predictions  of  the  electrical  energy  produced  by  the  PV 
array  are  based  on  a  GSI  prediction  model,  which  was  developed 
as  described  in  Section  3.2.  For  the  development  of  the  model,  we 
used  data  that  have  been  recorded  over  the  course  of  two  years  by 
the  meteorological  station  located  at  the  NTUA  campus,  Athens, 
Greece.  Fig.  7  shows  the  prediction  of  the  PVOut  three-day- 
profile  available  at  the  beginning  of  the  simulation.  The  PVOut 
prediction  corresponds  to  the  sequence  is  4  — >►  5  — >  4,  i.e.  partly 
cloudy  weather  conditions  for  the  first  day,  few  clouds  conditions 
for  the  second  day  and  partly  cloudy  conditions  for  the  third  day. 
The  equations  concerning  the  electrolyzer  and  the  fuel  cells  have 
been  provided  by  a  manufacturer  located  in  Greece  [27].  The 
equations  take  into  account  the  parasitic  loads  that  are  present 
due  to  the  auxiliary  subunits  in  each  particular  system. 

/ELOut (k)\ 

ELIn(k)  =  3.348  ELOut(k)  +  1.367  In  f  ^  ) 

+  3.845  x  10"5  ELOut(A:)2  (22) 


Fig.  6.  Reference  load  profile  for  October  reference. 
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Fig.  7.  Prediction  of  the  PV  output  three-day -pro file. 


Table  1 

Hybrid  system  specifications 


First  scenario 

Second  scenario 

PV  array 

4kWp 

4kWp 

Electrolyzer 

0.25  Nm3  h-1 

0.5  Nm3  h-1 

Fuel  cell  stack 

1.2  kW 

1.2  kW 

Metal  hydride  tanks 

5  Nm3 

5  Nm3 

FCOut(ifc)  =  -0.273  x  1 0-9  FCIn  (k  f 

+7.726  x  10“7  FCIn  (k)3  —  (23) 

0.916  x  10-3  FCIn  (k)2  +  1 .745  FCIn  (k) 

As  explained  in  Section  3.3  the  inverter  efficiency  is  assumed 
constant  and  equal  to  89.3%.  The  cost  for  purchasing  electrical 
energy  is  0.09  €  kWh-1  whereas  it  is  assumed  that  the  profit 
from  selling  electrical  energy  to  the  grid  is  0.05  €  kWh-1.  A 
lower  bound  of  1  Nm3  (20%  of  the  hydride  tanks  storage  capac¬ 
ity)  is  imposed  on  the  hydrogen  inventory  level.  Due  to  the  fact 
that  we  use  hydride  tanks  of  small  sizes  in  our  simulations,  the 
cost  for  removing  or  providing  heat  to  the  tanks  is  assumed  neg¬ 
ligible  and  therefore  that  parameters  CostCool  and  CostHeat  are 
set  equal  to  zero. 

Two  scenarios  regarding  the  specification  of  the  RESHS  were 
simulated.  The  particular  details  are  summarized  in  Table  1 
which  shows  the  production  capacities  of  the  three  main  pro¬ 
cesses  (PV  array,  electrolyzer,  fuel  cell)  and  the  storage  capacity 
of  the  metal  hydride  tanks.  The  only  difference  between  the  two 
scenarios  is  the  maximum  production  rate  of  the  electrolyzer, 
which  is  lower  in  the  first  scenario.  As  initial  and  final  values 


for  the  hydrogen  inventory  level,  50%  and  20%  of  the  storage 
capacity  were  assumed  in  both  scenarios. 

Fig.  8  presents  graphically  the  solutions  of  the  optimization 
problems  which  were  formulated  for  the  two  scenarios.  The  solu¬ 
tions  are  presented  in  the  form  of  optimal  profiles  (plots  (a)-(h)) 
of  all  input  and  output  variables  over  the  future  time  horizon 
of  72  h.  The  first  scenario  is  represented  by  green  lines  and  the 
second  scenario  by  blue  lines. 

Regarding  the  first  scenario,  we  notice  that  during  the  hours 
of  peak  generation  of  electrical  energy  from  the  PV  array,  the 
electrolyzer  operates  close  to  its  maximum  production  rate.  Dur¬ 
ing  the  same  period,  electricity  from  the  grid  is  zero  whereas  the 
surplus  is  sent  to  the  grid.  Furthermore,  a  gradual  filling  of  the 
tanks  with  hydrogen  is  observed  due  to  the  operation  of  the  elec¬ 
trolyzer.  This  inventory  build-up  is  used  for  producing  electricity 
in  the  fuel  cell  at  subsequent  periods  when  there  is  high  energy 
demand  but  low  or  zero  generation  of  electrical  energy  from  the 
PV  array.  So,  in  the  evening  hours,  when  energy  demand  is  high 
but  solar  energy  is  no  longer  available,  the  demand  is  mainly  ful¬ 
filled  by  the  fuel  cells,  while  small  amounts  of  electrical  energy 
are  purchased  from  the  grid.  Comparing  the  optimal  profiles  of 
the  two  scenarios,  it  is  clear  that  the  specification  is  critical  for 
the  performance  of  the  RESHS.  Due  to  the  increased  produc¬ 
tion  capacity  of  the  electrolyzer,  the  cost  for  purchasing  electrical 
energy  from  the  grid  is  lower  in  the  second  scenario.  Table  2  sum¬ 
marizes  the  results  for  the  first  scenario  regarding  the  amounts 
of  energy  that  are  purchased  from  the  grid  or  sold  throughout 
the  entire  time  horizon  and  the  corresponding  objective  function 
value. 

In  the  results  presented  so  far,  we  have  assumed  that  the  tem¬ 
poral  profiles  of  the  input  variables  (i.e.  PVOut  and  CS)  are 
deterministic,  i.e.  their  true  values  match  exactly  the  initial  pre¬ 
dictions  throughout  the  prediction  horizon.  Obviously,  perfect 
prediction  cannot  be  achieved  in  real-life  applications.  In  order 
to  examine  this  issue,  we  generated  modified  CS  temporal  pro¬ 
files  by  introducing  noise  to  the  available  predictions  whereas 
for  the  PVOut  we  used  true  profiles  as  they  have  been  recorded 
for  three  subsequent  days  in  October.  In  particular,  the  value  of 
CS (&+/),  at  each  future  time  instance  /=  1,. .  .,P  was  distorted 
by  adding  a  random  number,  chosen  from  a  uniform  distribution 
in  the  interval  [— 0. 1  -CS(k  +  /)  0. 1  -CS(k  +  /)] .  Regarding  the  true 
PVOut  values,  the  recorded  data  [24]  for  the  three  subsequent 
days  correspond  to  partly  cloudy  (state  4),  clear  sky  (state  6)  and 
heavy  clouds  (state  2)  conditions  respectively.  Fig.  9  presents 
the  PVOut  recorded  data  and  the  modified  CS  profile  together 
with  the  available  predictions.  It  is  clear  that  as  far  as  PVOut 
is  concerned,  we  introduce  systematic  error  to  the  available 
predictions,  which  increases  with  the  period  into  the  future. 


Table  2 

Purchased/sold  energy  and  objective-function  value  using  deterministic  and  updated  decision  strategies 


Total  energy  from 

Total  energy  to  the 

Objective 

the  grid  (kWh) 

grid  (kWh) 

function  (€) 

Deterministic  schedule  using  predicted  PVOut  and  CS  temporal  profiles 

9.7808 

4.9563 

0.63 

Deterministic  schedule  using  true  PVOut  and  CS  temporal  profiles 

27.8317 

17.2793 

1.64 

Updated  schedule  using  true  PVOut  and  CS  temporal  profiles 

13.2176 

7.9778 

0.79 
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Given  the  difference  between  true  and  predicted  PVOut  and 
CS  profiles,  an  obvious  question  emerges:  how  might  this  dis¬ 
crepancy  deteriorate  the  efficiency  of  the  decision  policy  that 
was  obtained  as  the  solution  to  the  deterministic  optimization 
problem?  In  order  to  answer  that  question,  we  used  again  the 
deterministic  schedule  that  was  determined  for  Scenario  1,  but 
now  considering  the  true  instead  of  the  simulated  PVOut  and  CS 
temporal  signals.  During  each  time  period  additional  shortages 
of  energy  were  satisfied  by  purchasing  it  from  the  grid,  while 


additional  surpluses  of  energy  were  sold  to  the  grid.  The  results 
are  shown  with  green  lines  in  Fig.  10.  The  amounts  of  energy  that 
are  purchased  from  the  grid  or  sold  throughout  the  entire  time 
horizon  are  shown  in  Table  2  and  compared  with  those  obtained 
when  zero  noise  was  assumed.  Clearly,  both  quantities  are  con¬ 
siderably  higher  when  the  true  temporal  PVOut  and  CS  profiles 
are  taken  into  consideration.  The  total  cost  (value  of  the  objec¬ 
tive  function)  increases  substantially  (160%).  This  observation 
actually  means  that  efficiency  becomes  lower  when  the  sched- 


Fig.  8.  Decision  strategy  for  two  different  scenarios  of  a  RESHS:  (a)  electricity  from  the  grid;  (b)  electrolyzer  consumption;  (c)  electrolyzer  production;  (d)  fuel  cell 
consumption;  (e)  fuel  cell  production;  (f)  electricity  to  the  grid;  (g)  direct  transferred  power  and  (h)  inventory  (scenario  1  and  scenario  2  are  represented  by  green 
and  blue  lines,  respectively).  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 
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Fig.  8.  ( Continued ). 


ule  obtained  from  the  solution  of  the  deterministic  optimization 
problem  is  applied  to  realistic  situations  containing  noisy  PVOut 
and  CS  temporal  profiles. 

In  the  rest  of  the  section,  we  will  explore  the  benefits  of  using 
the  rolling  horizon  concept  for  improving  the  decision  strategy 
that  has  been  determined  so  far.  The  rolling  horizon  philosophy 
is  based  on  the  fact  that  at  each  discrete  time  instance,  the  formu¬ 
lation  of  the  optimization  problem  can  be  updated  with  the  true 
current  state  of  the  system  (the  hydrogen  inventory  level)  and 
with  more  reliable  predictions  regarding  the  future  PVOut  and 
CS  temporal  profiles.  In  our  particular  application,  we  assume 
that  at  each  time  instance  we  know  the  true  values  of  PVOut 
and  CS  only  for  the  next  time  period,  while  for  the  rest  of  the 
future  horizon  we  keep  using  the  initial  PVOut  and  CS  predic¬ 
tions.  This  assumption  is  rather  conservative,  since  in  real  life 
situations  we  can  most  probably  improve  our  predictions  for 
the  entire  future  time  horizon.  However,  as  we  will  see  in  the 
sequel,  our  assumption  is  good  enough  for  demonstrating  the 
effectiveness  of  the  rolling  horizon  strategy. 

Following  the  previously  described  information  update  pro¬ 
cedure  at  the  end  of  each  discrete  time  period,  the  complete 


optimization  problem  is  reformulated  and  solved.  The  optimal 
control  strategy  resulting  from  the  rolling  horizon  philosophy 
is  shown  in  red  lines  in  Figure  10,  where  we  can  compare 
it  with  the  decision  strategy  produced  by  the  solution  to  the 
deterministic  problem  (shown  in  green  lines).  We  observe 
slight  modifications  in  the  electrolyzer  operation  schedule  for 
the  first  two  days  and  major  modifications  during  the  third 
day.  For  the  fuel  cell  system  considerable  modifications  are 
observed,  which  result  to  a  quite  different  inventory  temporal 
profile  and  to  reduced  amounts  of  purchased  and  sold  elec¬ 
tricity  (see  Table  2).  The  total  cost  is  quite  close  to  the  one 
calculated  by  using  the  predicted  and  not  the  actual  PVOut 
and  CS  temporal  profiles  (there  is  only  a  25%  increase).  It  is 
also  useful  to  compare  this  result  with  the  cost  for  covering  all 
the  consumer  needs  by  purchasing  electrical  energy  from  the 
grid,  which  is  3.16  €  (i.e.  only  25%  of  the  consumer  needs  is 
covered  by  purchasing  electrical  energy  from  the  grid).  There¬ 
fore,  the  adoption  of  the  rolling  horizon  concept  can  improve 
considerably  the  efficiency  of  the  hybrid  system  in  real-life  sit¬ 
uations,  where  PVOut  and  CS  predictions  obviously  contain 
errors. 


Fig.  9.  True  and  predicted  temporal  profiles  for  (a)  PVOut  and  (b)  CS. 
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Fig.  10.  Deterministic  and  updated  decision  strategies  for  scenario  1  with  true  PVOut  and  CS  profiles  (deterministic  and  updated  decision  strategies  are  represented 
by  red  and  green  lines  respectively).  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 


5.  Conclusions 

In  this  work,  a  complete  framework  for  managing  renew¬ 
able  energy  systems  with  hydrogen  storage  is  proposed.  The 
framework  is  built  in  several  steps  that  may  be  summarized  as 
follows: 

•  Development  of  a  prediction  model  for  the  GSI  daily  pro¬ 
file  on  horizontal  surfaces  based  on  the  RBF  neural  network 
architecture; 

•  Forecasting  GSI  on  the  tilted  surface  of  the  PV  array,  using 
the  following  information:  the  GSI  prediction  on  horizontal 
surfaces,  astronomical  and  geographical  data  and  the  slope  of 
the  PV  array; 

•  Estimating  the  electrical  energy  produced  by  the  PV  array, 
given  the  ambient  temperature  and  the  GSI  on  the  tilted  sur¬ 
face; 

•  Development  of  a  model  that  realistically  describes  the  per¬ 
formance  and  the  constraints  of  the  hybrid  system; 

•  Formulation  and  solution  of  an  on-line  optimization  problem 
that  is  used  as  a  decision  making  tool  regarding  the  operation 
of  the  system.  The  formulation  takes  into  account  updated 
estimations  of  the  photovoltaic  power  generation  over  a  future 
prediction  horizon  and  a  profile  of  the  energy  demand  over 
the  same  time  horizon. 

In  conclusion,  the  proposed  model  may  prove  to  be  a  very 
useful  tool  for  optimal  decision  making  in  hybrid  power  gener¬ 
ation  systems  which  combine  RES  and  hydrogen  technologies. 
It  should  also  be  mentioned  that  the  actual  cost  of  all  individual 
components  of  such  a  hybrid  system  is  continuously  decreasing, 


since  renewable  energy  technologies  become  mature.  Hybrid 
renewable  energy  systems  are  expected  to  be  competitive  to  the 
conventional  power  generation  systems  in  the  near  future  and 
thus  optimization  of  their  operation  is  of  particular  interest. 
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